#include <limits>
#include "point.h"
#include "adaptive.h"

const std::size_t size = std::numeric_limits<double>::digits / 2 + std::numeric_limits<double>::digits % 2;

struct pair_double
{
   double low, high;

   pair_double(double a = 0, double b = 0): low(a), high(b)
   {
   }
};

namespace
{
   pair_double sum(double a, double b)
   {
      double low = a + b;

      double bv = low - a;
      double av = low - bv;

      double br = b - bv;
      double ar = a - av;

      double high = ar + br;

      return pair_double(low, high);
   }

   template <std::size_t s>
   pair_double split(double a)
   {
      double c = ((1LL << s) + 1LL) * a;
      double ab = c - a;

      double high = c - ab;
      double low = a - high;

      return pair_double(low, high);
   }

   pair_double mul(double a, double b)
   {
      double low = a * b;

      pair_double a_split = split<size>(a);
      pair_double b_split = split<size>(b);

      double e1 = low - a_split.high * b_split.high;
      double e2 = e1 - a_split.low * b_split.high;
      double e3 = e2 - b_split.low * a_split.high;

      double high = a_split.low * b_split.low - e3;

      return pair_double(low, high);
   }
   
   template <std::size_t m>
   void grow_expansion(double* e, double b, double* h)
   {
      double q = b;

      for(std::size_t i = 0; i < m; i++)
      {
         pair_double temp = sum(e[i], q);
         h[i] = temp.high;
         q = temp.low;
      }

      h[m] = q;
   }

   template <std::size_t m, std::size_t n>
   void expansion_sum(double* e, double* f)
   {
      for(std::size_t i = 0; i < n; i++)
         grow_expansion<m>(e + i, f[i], e + i);
   }
   
   template <std::size_t n>
   int adaptive_turn_sign(double * e)
   {
      for(int i = n - 1; i >= 0; i--)
      {
         if(e[i] > 0)
            return 1;

         if(e[i] < 0)
            return -1;
      }

      return 0;
   }
};

int adaptive_turn(const point& a, const point& b, const point& c)
{
   double mult[12];
   pair_double temp;

   temp = mul(b.x, c.y);
   mult[0] = temp.high;
   mult[1] = temp.low;

   temp = mul(-b.x, a.y);
   mult[2] = temp.high;
   mult[3] = temp.low;     

   temp = mul(-a.x, c.y);
   mult[4] = temp.high;
   mult[5] = temp.low;

   temp = mul(-b.y, c.x);
   mult[6] = temp.high;
   mult[7] = temp.low;

   temp = mul(b.y, a.x);
   mult[8] = temp.high;
   mult[9] = temp.low;

   temp = mul(a.y, c.x);
   mult[10] = temp.high;
   mult[11] = temp.low;
      
   expansion_sum<2, 2>(mult, mult + 2);
   expansion_sum<2, 2>(mult + 4, mult + 6);
   expansion_sum<2, 2>(mult + 8, mult + 10);
   expansion_sum<4, 4>(mult, mult + 4);
   expansion_sum<8, 4>(mult, mult + 8);      

   return adaptive_turn_sign<12>(mult);
}